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The dynamics of the collective excitations of a lattice Bose gas at zero temperature is systemati- 
cally investigated using the time-dependent Gutzwiller mean-field approach. The excitation modes 
are determined within the framework of the linear-response theory as solutions of the generalized 
Bogoliubov-de Gennes equations valid in the superfluid and Mott-insulator phases at arbitrary val- 
ues of parameters. The expression for the sound velocity derived in this approach coincides with 
the hydrodynamic relation. We calculate the transition amplitudes for the excitations in the Bragg 
scattering process and show that the higher excitation modes give significant contributions. We 
simulate the dynamics of the density perturbations and show that their propagation velocity in 
the limit of week perturbation is satisfactorily described by the predictions of the linear-response 
analysis. 
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I. INTRODUCTION 

Studies of excitations of ultracold atoms in optical lat- 
tices play an important role for understanding of their 
physical properties and dynamical behavior. For in- 
stance, celebrated quantum phase transition from the 
superfluid (SF) into the Mott-insulator (MI) [l[ is ac- 
companied by the opening of the gap in the excitation 
spectrum In the deep optical lattices the system of 
interacting bosons is satisfactorily described by the Bosc- 
Hubbard model [H, 0- Since the model is not integrable, 
exact analytical results can be obtained only in few spe- 
cial cases like in the limit of weakly interacting gas Q or 
for hard-core bosons in one dimension Q. However, in 
general exact results can be obtained only with the aid 
of numerical methods. 

Exact numerical results for the spectrum of low-energy 
excitations were obtained by means of diagonalizations of 
the Bosc-Hubbard Hamiltonian 

Although numer- 
ical diagonalization can be really done only for rather 
small lattices which are far from the thermodynamic 
limit, this allows to capture all the main characteristic 
features of realistic systems. Numerical results for larger 
systems of the same size as in real experiments with ul- 
tracold atoms [8j|_have been obtained by quantum Monte 
Carlo methods J9| which allow to compute the spectral 
properties [lol lll| . Simulation of the real-time dynam- 
ics of quantum systems within quantum Monte Carlo is 
also possible (see, e.g., [l2j) but has not been yet per- 
formed for Bose systems. Ground-state properties as 
well as the real-time dynamics of one-dimensional sys- 
tems subject to external perturbation were studied, e.g., 
in Refs. [13T - [lq ] by the powerful numerical density-matrix 
renormalization group (DMRG) method giving an access 
to the excitation spectrum [17| but so far restricted to 
one-dimensional systems. 

Mean-field theories in the dimensions higher than one 



allow self-consistent study of the excitations and dynam- 
ics in the lowest order with respect to quantum fluctu- 
ations. In a weakly interacting regime, the atoms are 
fully condensed and the system is satisfactorily described 
by the time-dependent discrete Gross-Pitacvskii equation 
(DGPE). The excitation spectrum can be calculated us- 
ing Bogoliubov-de Gennes (BdG) equations [l8l - [20l |. It 
has a form of the Goldstone mode characterized at large 
wavelengths by the sound velocity c. The latter is related 
to the compressibility k, the effective mass m* and the 
condensate density |i/>| 2 through the relation 



c = \/\ip\ 2 /nm* 



(1) 



The same expression for the sound velocity can be also 
derived from the Bogoliubov theory in the operator for- 
malism 0] and from the hydrodynamic approach [l|, l2l| . 

In the strongly interacting regime, condensate fraction 
becomes suppressed and for commensurate fillings the 
system can undergo a transition from the SF into the MI 
state 0, Hi], [22[. In this regime, DGPE is not valid and 
on the mean-field level must be replaced by more general 
Gutzwiller equations (GE) which are exact for a gas of 
infinite dimensions |l|, l23l425l | . GE were successfully used 
to study various phenomena like creation of molecular 
condensate by dynamically melting a MI (26j , dynamical 
transition from the SF to Bose-glass phase due to con- 
trolled growing of the disorder 1271. t he gas dynamics in 
time-dependent lattice potentials |28| . transport of cold 
atoms induced by the shift of the underlying harmonic 
potential 1291 , dynamics of metastable states of dipolar 
bosons |30l|. and the soliton propagation [3l| . 

Excitations above the ground state described by the 
GE were studied using the random phase approximation 
(RPA) [32l436| , the Schwinger-boson approach [37j , time- 
dependent variational principle with subsequent quan- 
tization 13811 . Hubbard-Stratonovich transformation [39L 
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representation of the Bose-Hubbard 
standard-basis operator method [13, EH, and 
Ginzburg-Landau theory [44|. All these methods show 
that in the MI phase the spectrum of excitations consists 
of the particle- and hole-modes, both with nonvanish- 
ing gaps [13 EH • In the SF phase, the lowest branch is 
a gapless Goldstone mo de [37l HH, ES EH , while higher 
branches are gapped [37l |46|. However, different approx- 
imations used in the calculations by different methods do 
not always lead to the same final results. For instance, in 
Ref. it was checked numerically that the RPA gives 
the same result for the sound velocity as Eq. AT}, whereas 
the analytical expressions derived in Refs. [4fJ, |47j differ 
from that. 

Excitations above the Gutzwiller ground state can be 
also investigated using generalization of the BdG equa- 
tions directly derived from the GE within the framework 
of the linear response theory. This method was used for 
the lattice Bose gas with short-range 0, E3] as well as 
long-range [30l |4o| interactions. This approach, which 
was not so far widely used to study the lattice Bose gas 
allows to obtain results consistent with other mean-field 
approaches mentioned above and to study the ground 
state, stationary excitation modes as well as the dynam- 
ics of the gas on equal footing. 

The excitations can be probed in experiments on in- 
elastic light scattering (Bragg spectroscopy) which pro- 
vide an information on the dynamic structure factor 
and one-particle spectral function. Recently such ex- 
periments were carried out with ultra-cold rubidium 
atoms in o ptic al lattices of different dimensions in the 
SF phas e El El, EH and across the SF-MI transi- 
tion l5l|. Theoretical analysis has been developed 
using exact diagonalization of the Bose-Hubbard Hamil- 
tonian [52l l53j, q uantum Monte Carlo simulations in 
one dimension [1(1 53j. perturbation theory valid deep 
in the MI phase [52j. hvdrodvnamic theory exact 
Bethe-ansatz solution of the Lieb-Liniger model [55| , ex- 
tended fermionization (53|, RPA (HJ- Analogous stud- 
ies were also performed within the Gutzwiller approxi- 
mation [HI, H3EH- However, previous calculations are 
valid only for the MI [H [13 or close to it [13 • It was 
argued that the second excitation branch in the SF can- 
not be detected by the Bragg spectroscopy in the linear 
regime [38| . In Ref. [46|, the calculations beyond the 
linear response theory taking into account the harmonic 
trapping potential and the finite duration of the probing 
Bragg pulse were performed which are in good quanti- 
tative agreement with the experimental data reporting 
the observation of the second excitation branch. In the 
present work, we will study in details the possibility to 
observe the second excitation mode in a homogeneous 
lattice in the linear-response regime. 

Experimentally sound waves can be observed with the 
aid of an external potential which creates a density per- 
turbation of the gas. Corresponding numerical simula- 
tions for the lattice gas were performed on the basis of 
the DGPE [1H and for soft-core bosons in ID making 



use of the DMRG method 14| . The sound velocity ex- 
tracted from these simulations is in perfect agreement 
with Eq. ([1]) in the case of the DGPE [2(| and has a cor- 
rect asymptotic behavior in ID in the limits of weak and 
strong interactions 1J] , where analytical expressions are 
known. 

The purpose of this paper is to give a comprehen- 
sive self-consistent description of the collective excita- 
tions as well as experimental techniques for their observa- 
tion within the time-dependent Gutzwiller ansatz which 
is gapless and satisfies the basic conservation laws, in 
particular, f-sum rule. Solution of GE allows not only 
to obtain the excitation dispersion relations but also to 
calculate the transition amplitudes in the Bragg scatter- 
ing process. We present a derivation of Eq. ([1]) from 
the generalized BdG equations. Furthermore, the time- 
dependent approach is also used to investigate the sound 
wave propagation in the case of a stronger perturba- 
tion generated by switching off a local potential. In 
this way, we can determine the speed at which this per- 
turbation propagates and compare with the predictions 
of the linear-response theory. We emphasize that the 
Gutzwiller ansatz is the only approximation used in the 
present work and the results are valid in the whole range 
of parameters. 

The paper is organized in the following manner. In 
SecHH we present the time-dependent GE. Their ground 
state solutions are discussed in Sec. IIIIl In Sec. IIV1 we 
determine the spectrum of collective excitations using 
Gutzwillcr-Bogoliubov-de Gennes equations. Sec. [V] is 
devoted to the Bragg scattering. In Sec. I VII we simu- 
late the sound-wave propagation. The conclusions are 
presented in Sec. I VIII 



II. THE TIME-DEPENDENT GUTZWILLER 
ANSATZ 

We consider a system of ultracold interacting bosons in 
a e?-dimcnsional isotropic lattice described by the Bose- 
Hubbard Hamiltonian 



it; 



H = -J^2Y1 { a i ai 

1 a=l 



h.c. 

it 



l 



(2) 



where the site index 1 is a d-dimensional vector, e a is a 
unit vector on the lattice in the direction a, J is the tun- 
neling matrix element, U is the repulsive on-site atom- 
atom interaction energy, and /i the chemical potential. 
The annihilation and creation operators at site 1, ai and 
af, obey the bosonic commutation relations. The mo- 
mentum operator 



a ^2 ("l 



h.C. 
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where Pq is a constant determined by the parameters 
of the periodic potential creating the lattice, does not 
commute with the Hamiltonian ([2]) due to the interaction 
term. Instead, the quasi-momentum operator defined as 

p = ]T ka[a k , a k = £ e^%/L d / 2 , 

k 1 

where k a = 2nn a /L, n a = 0, ...,L — 1, with L being 
the number of lattice sites in each spatial direction, com- 
mutes with the latter. 

Our analysis employs the Gutzwillcr ansatz. Thereby, 
eigenstates of the Hamiltonian @ are taken as tensor 
products of local states 



(3) 



where |n)i is the Fock state with n atoms at site 1. Nor- 
malization of the |si) imposes 

oo 

EM 2 = i- 

The mean number of condensed atoms in this model is 
given by \ipi\ , where 



(4) 



is the condensate order parameter. One can easily show 
that l^i | 2 cannot be larger than the mean occupation 
number 



quasi-momentum, total energy, and the total number of 
particles remain constant in time. 

As it follows from the form of the state , Gutzwiller 
approximation neglects quantum correlations between 
different lattice sites but takes into account on-site quan- 
tum fluctuations. This appears to be enough for satisfac- 
tory description of the SF-MI quantum phase transition. 
Due to the fact that in the equations of motion (j6]) the co- 
efficients cin for different sites are coupled to each other, 
Gutzwillcr ansatz can be also used to study the dynamics 
of excitations. 

^From ((6J, we deduce the following equation for the 
order parameter: 



i— — 

dt 



a=l 
oo 

+ Uj2( n - 1 )V^cl n _ 1 c hn . (7) 

n=Q 

This equation becomes closed if we assume the coherent 
state 



(8) 



C, n = dff = cxp H^| 2 /2) W/y/ni 



which is an exact solution of Eq. ^ for U = 0. In this 
case, the replacement of the last term of Eq. ([7]) by this 
distribution leads to the result 

oo 

so that we recover the DGPE valid for small U/J. 



^2n\c ln \ 2 . 

n=l 

Minimization of the functional 

oo 

ih Y2(ci n d t Cln ~ Cl n d t C* ln ) - (H) 



n=0 



leads to the system of GE [23, |56|: 
. dc\ n 



ih- 



dt 



Hi " c lri , 



— n(n — 1) — [in 



5 n '. 



jVri/6n>,n+l (i'l+ ea + V'l-eJ 
ct=l 
d 

J \fn5 n>n i+i ^ (in+ea, + ik-e a ) 



(5) 



(6) 



The Gutzwiller approximation is conserving since these 
equations do not violate conservation laws of the origi- 
nal Bosc-Hubbard model. The expectation values of the 



III. GROUND STATE 

In the ground state, the coefficients c\. n do not depend 
on the site index 1 so that the solution has the form 



ci,„(f) = <4 0) exp (—iw t) , 



(9) 



so that (n\) = (n), where n = J2i n i/L d - The coeffi- 
cients Cn are calculated numerically solving Eq. ([5]) by 
means of exact diagonalization in the same manner as 
in Refs. [H, |3{|. The results are shown in Figs. [T] and 

[5] has a broad distribution in the SF phase, where 
the corresponding tfj^ = ipi defined by Eq. Q does not 
vanish. In the MI phase, however, 



(10) 



resulting in = 0. ujo is determined in both phases 
from: 



hu a = -4dJV' (0)2 + J2 



— n(n — 1) — fin 



4 0)2 - (11) 
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FIG. 1. Ground-state solutions for the atomic distribution 
c4 ■ The scaled chemical potential fi/U = 1.2 and the tun- 
neling rates 2dJ/U: 0.7 (i), 0.5 (ii), 0.3 (iii), 0.15 (iv), and 
0.05 (v). The lines connecting the dots are to guide the eye. 
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FIG. 2. Comparison of cL 0) (1,2) with c c ° h (a,b). The scaled 
chemical potential fi/U = 1.2 and the tunneling rates 2dJ/U: 
0.5 (l,a), 50 (2,b). 



In Figs. [2] and [3j we compare the coefficients cl? ob- 
tained by numerical solution of Eqs. ((BJ) with the corre- 
sponding results for the coherent state ((5]). The value 
of V (0) for Eq. © was calculated according to Eq. (HJ) 

using ci° . These coefficients converge towards a coher- 
ent state distribution for increasing J/U according to a 
power law (Fig. [3]), thus, justifying the use of the DGPE 
in this limit. 

In the numerical calculations presented in this section 
and later on, n was restricted by some finite N (c„ = 
for n > N). The cut-off number of atoms N was chosen 
large enough such that its influence on the eigenstates is 
negligible. For example, for the plots shown in Fig. [1] it 
was enough to use N = 10, but for Figs. EJ MN = 500 
was used. 



IV. EXCITATIONS 



We consider small perturbation of the ground state 
c\ n {t) = ci? + q„ (t) + ... cxp (-iu> a t), where 



FIG. 3. Deviation of the exact ground state from the coherent 
state for the scaled chemical potential fi/U = 1.2. 



Substituting this expression into GE and keeping only 
linear terms with respect to itkn and «kn, we obtain the 
system of linear equations [30| 



frcui. 



"k 
' 7 k 



"k 
' 7 k 



(13) 



where and arc infinite-dimensional vectors with 
the components Ukn and v^n (n = 0, 1, . . . ), respectively. 
Matrix elements of and £?k have the form 



+ — n(n — 1) — fin — Hloq S n ', n 



+ Jn V n' c. 



n-l Si'-l 



j->nn 



+ Jn Vn' + l c, 



(o) jo) 

n-l c n'+l 



where Jk = 2d J — ek with 



£k 



4 J ^ sii 



a=l 



(14) 



being the energy of a free particle. This system is valid 
for both phases and generalizes the BdG equations pre- 
viously derived for coherent states. The dependence on 
vector k is determined by the variable 



1/2 



(15) 



which varies from to 1. For small |k|, x « |k|/(2\/d). 

The energy increase due to the perturbation is given 
by |H 



(12) 



AE = hwk (\u k \ 2 - \v k \ 



(16) 
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Formally, Eqs. ([131 have solutions with positive and 
negative energies iftcjk, which are equivalent because 
Eqs. (|12| . (|16|) arc invariant under the transformation 
cjk — > — Wk, k — > k, Mk — > «k — > ""k, so that only 
solutions with the positive energies will be considered in 
the following. The eigenvectors are chosen to follow the 
orthonormality relations 



l k,A' 



Perturbation (fl2l) creates plane waves of the order pa- 
rameter ip\{t) = ip^ + t/ij (t), where 



Wk = ^ (c^Ukn + C^Vk.ri-l) , 



(17) 



c4 0) Uk,„-i 



(0) 



and 



The perturbations for the total density and the con- 
densate density are given by 

(ni)(t) = (h) + [^ke*^ 1 -'**) + ex.] , (18) 

Ak = 2J <4 0) 7J (Ufcn + «kn) , 

n 

|^(i)| a = V (0)2 + [Ske^' 1 -^ + c.c] , (19) 

n 

+ C^ 0) (Uk,n-1 + Uk,n-l) 

In what follows we consider the properties of the excita- 
tions in the MI and SF phases. Although the results for 
the MI are not new, we would like to present those for 
completeness. 

A. Mott insulator 

For the MI phase, the coefficients ch^ have a simple 
analytical form (fT0|). The eigenvalue problem for the 
infinite-dimensional matrices (fT3"|) reduces to the diag- 
onalization of two 2 x 2-matriccs which couple Uk,n -i to 
Vk,n +i an d Uk,n +i to i>k,n -i> respectively. The lowest- 
energy excitation spectrum consists of two branches with 
the energies 



lip 



± 



[/no 



2 



(20) 



The same result was obtained using Hubbard- 
Stratonovich transformation [39| and within the 
Schwingcr-boson approach [37| . 



These two branches are shown in Fig. [4] and display 
a gap. According to Eqs. p^|) . (fTS)) , no density wave is 
created in the two modes, although the order parame- 
ter does not vanish [see Eq. ([T7|) ]. Eq. can be also 
rewritten in the form 

?kj k + = £k P - Li , 
?kj k - = -e kh + fi . 

Therefore, the solutions with index and '— ' have 
the meaning of the particle and hole excitations, respec- 
tively H3- 

Other solutions of Eq. (fi"3]) are independent of k with 
the energies 

hu x = [A (A - 1) - n (n - 1)] - n(X - n ) , (21) 

They are denoted by A which are non-negative integers 
different from no,no i 1- If no is the smallest integer 
greater than n/U , the excitation energies are always pos- 
itive. The eigenvectors of these modes have the form 
Uk n x — 5n.\, Vknx = 0, and the amplitudes of all the 
waves defined by Eqs. ([TT|). (p~8|) . ([19]) vanish. 

The boundary between the SF and MI phases is deter- 
mined from the disappearance of the gap in the excitation 
spectrum, i.e., when wo- = 0. Under this condition, we 
recover the critical ratio H : 



2d(J/U) c 



(n - li/U){li/U -no + 1) 



1 + n/U 
It takes its maximal value when 

2d(j/u)™ ax = {v^+T-v^y 

for a chemical potential given by 



{H/U) c = v^ioN + l) - 1 



(22) 



(23) 



(24) 



For J/U > (J/U) c , the lowest frequency u; _ in Eq. (|2"0]) 
becomes negative leading to a negative expression for 
Eq. (Til)]) , so that the Mott-phase solution (TITJ)) does not 
correspond to the ground state anymore. 

The excitation spectrum has interesting features on 
the boundary between the MI and SF. For (J/U) c = 
(J/U)™ ax , the excitation energies (j2U)) can be rewritten 
as 



k± 



-,1/2 



y/no(no + l)^k + -f 



± 



fk 



(25) 



For small |k|, the two branches are degenerate and have 
linear dependence Wk± = c* ip |k| with the sound velocity 



U 



y/(J/U)™ x [no(no + l)] 



i/4 



(26) 



expressed in the units of number of sites per second. For 
other points on the boundary, i.e., {J/U) c < (J/U)™ ax , 
no degeneracy appears and the sound velocity vanishes 
leading to the quadratic dispersion Wk± ~ k 2 for small 
Ikl. 
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FIG. 4. First three branches Tilu^- (1), 7kdk+ (2) and 
Hlj\ = o (3) of the excitations spectrum of the MI for fj,/U = 1.2 
and 2dJ/U = 0.05, which corresponds to no = 2. 
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FIG. 6. Band structure of the excitation spectrum for (n) — 
0.5 (a), 1 (b). Shaded regions which are extremely narrow 
for higher bands show allowed excitation energies. Lower and 
upper boundaries of the bands correspond to x — and x = 1, 
respectively. 



FIG. 5. First three branches hto^^x (A = 1,2,3) of the exci- 
tations spectrum of the SF for fj,/U = 1.2 and 2dJ/U = 0.15. 
The straight dashed line represents the linear approximation 
with the sound velocity (|28l) . 



B. Superfluid 

In the SF phase, the eigenvalue problem (fl"3]) is solved 
using the numerical values of c„ for each J/U and [i/U. 
The energies of the lowest-energy excitations are shown 
in Fig. [5J The excitation spectrum consists of several 
branches which form a band structure shown in Figs. [51 [7] 
In contrast to the MI, the lowest branch has no gap. It is 
a Goldstone mode which appears due to the spontaneous 
breaking of the phase symmetry and, therefore, can be 
called a phase mode [38[ • As it is shown in Fig. [HJ the 
amplitude of the total-density wave is larger than the am- 
plitude of the condensate-density wave for this mode. A 
value greater than unity for the ratio Au/Bk means that 
the condensed part and normal part oscillate in phase. 

Higher modes (A > 2) do not exist in the formalism 
based on the DGPE. They have gaps A> = hu)o\ which 
grow with the increase of J (see Fig. [5]). As it is shown 
in Appendix [C] for large J they have an asymptotic form 



A x = 2dJX+ -A((n) + A 



1) 



(27) 



We note that only the first two lowest-energy branches 
have a strong dependence on k. For the second mode 



(A = 2), the amplitude of the total-density wave is 
much less than that of the condensate-density wave (see 
Fig. [5]) which means that the oscillations of the conden- 
sate and normal components are out-of-phase. How- 
ever, this does not necessarily mean that there is an 
exchange of particles between the condensate and nor- 
mal component J37l . |46| . Due to the reasons explained 
in Refs. [37], HH |4(f this type of excitations is called an 
amplitude mode. In the rest part of this section, we will 
study in more details the properties of the Goldstone 
mode. 

As shown in Appendix [XI the lowest-energy branch 



has a linear form Wk,i 
velocity given by 



Ikl for small k with the sound 



^(0) 



/h, 



(28) 



where n = g^- is the compressibility. This result proves 
that the Gutzwillcr approximation is gapless and coin- 
cides with Eq. (JJ). 

Fig. [§] shows the dependence of the sound velocity on 
fi and J calculated numerically using Eq. ([28)) . If we ap- 
proach the boundary of the MI, the sound velocity goes 
to zero everywhere except the tips of the lobes, where 
it is perfectly described by Eq. ([2l)]) . This behavior can 
be understood considering the properties of ip^ & n d k. 
If we approach the SF-MI transition from the SF part of 
the phase diagram, the order parameter tends always 
continuously to zero. The compressibility k reaches a fi- 




X 



FIG. 7. Band structure of the excitation spectrum for 

2dJ/U = 0.2 (a), 0.3 (b) versus the density. Shaded regions FIG. 8. A k /B k for fi/U = 1.2 and 2d,J/U = 0.15 (a), 1 (b). 
show allowed excitation energies. Lower and upper bound- 
aries of the bands correspond to x = and x = 1, respec- 



nite value at every point of the boundary except the tips 
of the Ml-lobcs where it tends continuously to zero such 
that the ratio 4> is finite. Therefore, the sound ve- 

locity vanishes at any point of the phase boundary except 
the tips of the lobes 0, HU . 

For a weakly interacting gas (U <C J), | -0 C°) | 2 « (n) 
and k ~ 1/U. In this limit, we recover the Bogoliubov's 
dispersion relation (see Appen dix [Bjl and the expression 
for the sound velocity [2fJ, [59[ 

cf = ^j2JU{h)/h . (29) 

The comparison of the exact numerical values of the 
sound velocity calculated according to Eq. (|28[) with the 
approximation (f2"9")l is shown in Fig.[TUJ As it is expected, 
the agreement is good at large J but for small tunneling 
rates the behavior of cf is completely different. 

In the opposite limit ( J <C U), the superfluid regions 
with the atomic densities no < (h) < no + 1 are confined 
in the regions /i_ < fi < where a_ is the upper 
boundary of the MI with n and u + is the lower boundary 
of the MI with no + 1. For J — > 0, fj, is a linear function 
of the density, i.e., 

M = M- + (M+ - M-) ((«■) - ™o) ■ (30) 
Using Eq. ([22]) up to the first order in J, we obtain 

/Lt± = Uno ± 2rfJ(n + 1) (31) 
from which we deduce 

k = l/[4dJ(n + 1)] ■ (32) 




FIG. 9. Sound velocity calculated numerically from Eq. (|28|) . 
Note the discontinuities at the points [(J /U)™ ax , (fJ,/U) c ] de- 
scribed by Eq. ([26)) . 

Using dE/d(h) = u and the condition that the energy 
at the phase boundaries E± = Un (n ± l)/2, the total 
energy becomes in this limit 

E = rjn ((n)-^±1) (33) 

- 2d J (n + 1) ((n) - n ) (no + l-(n}) . 

In order to deduce the order parameter, we use the rela- 
tion dE/dJ = -2d\ipW\ 2 . From Eq. we obtain 

^(°)| 2 = (n + l)«n) - no)(no + 1 - (n)) . (34) 

Eqs. (|3^|) . ([3^1) coincide with the ones obtained in Ref. [24[ 
but using a perturbation approach. Substituting (|32[) and 
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FIG. 10. Comparison of the sound velocity calculated numer- 
ically from Eq. ()28p (solid line) with the analytical expression 
(pi))! (dashed line) for (n) = 0.5. 



34)) into (f28|) . we obtain finally: 



c° = 2J(n + l)\/2d((n) - n Q )(n Q + 1 - (h))/h .(35) 

The sound velocity in the limit of small J vanishes at 
(n) = no, no + 1 and takes maximal values at (n) = uq + 
1/2. This qualitative behavior is the same as in the case 
of hard-core bosons in ID, where the sound velocity is 
given by c^ c = 2Jsin(7r(n))/?i for < (n) < 1 (see, 
e.g., 1). 

Close to the phase boundary, the sound velocity can 
be calculated analytically according to Eq. (f28]l within 
the fourth-order perturbation theory [HI]. The result 
for c s appears to be very long and cannot be displayed 
here. However, it describes all the qualitative features 
discussed above and reproduces Eq. ([26)l. 



V. TRANSITION PROCESSES 
A. Bragg scattering 

The Bragg scattering is a common experimental 
method to measure the excitation spectrum of an ultra- 
cold gas [18j . This process is induced by the following 
perturbation term in the Hamiltonian: 

H'(t) = ^ Vk, w cos(k • 1 - a;t)ajai , (36) 



where co and k are the frequency and the wavevector of 
the excitation, respectively. This perturbation changes 
the gas density according to: 

5m(t) = (Ai(f)> - (ft) = ^Sp^e^-^+c.c. (37) 

Up to the first order, this change is linear in the Bragg 
potential, i.e., 



where x(k, w) is the susceptibility. Within the Gutzwillcr 
approximation, we obtain 



x(k^) 



(39) 



'o 

r (0) 



A^ — huj Bk 
£>k A^ + Hlu 



r (o) 
'o 

r (0) 



where to = uj + iO and the vector components are defined 
as u °n = ncj? . After matrix inversion, the susceptibility 
can be rewritten in the more simple form 

2 \ ^ Xk.A^k.A 



ioy 



(40) 



k.A 



where A denotes various excitation branches associated 
to the eigenvalues u;k,A an d 



Xk : A 



? (0) 



(wk,. 



«k.A, 



(41) 



is the amplitude for the Bragg scattering of transition 
frequency Wk,A- As we see, Xk,A is nothing but the square 
of amplitude of the density wave defined by Eq. (fl8)) . 

For long wavelength, only the lowest mode is domi- 
nant since Xk,A vanishes for the higher modes. This is a 
consequence of the orthogonality between the eigenvec- 
tor components of the other mode and (u °\ — u °' 1 ) in 
the long- wavelength limit (see Appendix [5}. Since A = 1 
denotes the sound branch, we obtain the approximate 
expression 



k^O 



Xk.i^k.i 



h (oj + iO) 2 - 
The comparison of (|4"2")l with the identity 
K = -x(k = 0,^ = 0) , 



(42) 



(43) 



which follows from Eq. 
long wavelength 



, allows to deduce that for 



k->0 ft on I 

Xk,i = 



(44) 



<5pk,^ = x(k,w)Vk,u, , 



(38) 



The dependences of Xk,A on the variable x defined by 
Eq. (JT5")) for the excitation branches with A = 1,2,3 are 
shown in Fig. [TTJ For the chosen values of parameters, 
only two lowest branches display noticeable amplitudes. 
Similar results have been also obtained in Ref. [331 ■ How- 
ever, the calculations in Ref. (3?} are valid only close to 
the boundaries MI-SF because the occupation numbers 
n in Eq. ([3]) were restricted to n — hq, tiq ± 1. 

In Fig. [TJ] instead, we see that the amplitude for the 
third excitation branch as well as for the second one can 
become significant at certain densities. We would like to 
note that the f-sum rule is automatically fulfilled in our 
approach (see Sect. IV C")) in contrast to Refs. (33l. l37j. 

In the MI phase, the Gutzwiller approximation does 
not allow to observe any branches since Xk.A = 0. No 
Bragg response is possible, although the excitations ex- 
ist in the mean-field approach. In order to allow non- 
vanishing response, correlations between different sites 
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FIG. 11. Transition amplitudes Xk,A associated the transition 
frequency Wk,A for the lowest excitation branches (A = 1, 2, 3) 
and for /i/U = 1.2 and 2dJ/U — 0.15. The dashed line corre- 
sponds to the approximation ((44)) . 




0.6 




FIG. 12. (color online) Transition amplitudes Xk,A versus the 
density for the first excitation branches (A = 1,2,3) and for 
x = 1 and 2dJ/U = 0.2 (a), 0.3 (b). Dashed lines show the 
static structure factor S(k). 



should be included which goes beyond the Gutzwillcr 
approximation In such a description, excita- 

tions in the Bragg process are created as particle-hole 
pairs [H, [ID, H3|- However as pointed out in 25 1, this 
last process appears to be of the second order in the in- 
verse of coordination number z = 2d and, therefore, is not 
taken into account by the standard Gutzwillcr ansatz. 



B. One- particle Green's function 

The one-particle Green's function can be also deter- 
mined in the context of the Gutzwillcr approximation 
through the interaction term 



(45) 



which explicitly breaks the U(l) symmetry. This inter- 
action term induces a deviation in the order parameter 



ipt - (v> (0) r 



Hi 

HI 



(46) 



The proportionality term is the one-particle 2x2 matrix 
Green's function with the general expression 



£(k,a,) = £ 



uj + iO - Wk,A 



(47) 



where we define the matrix transition amplitude as: 

fj. , = & k ,A-&L (48) 



=k,A 



and 



ik,A 




lWk,n+l,A + y/nVk t n-l,\) 
ri+l,A + \AlUk,n-l,A) 



(0) 
Cn 



(49) 

Here, A denotes branches with both positive and negative 
energies. In the SF phase, the existence of the order pa- 
rameter hybridizes the one- and the two-particle Green's 
functions so that their poles are identical. On the other 
hand, in the MI phase, we note that the transitions for- 
bidden in the Bragg scattering become allowed in the 
interaction term Eq. (|45p . Indeed, the time-dependent 
Gutzwillcr approach allows to recover the results previ- 
ously established in the context of quantum field the- 
ory HU: 



9k,± 



where 



1 , (2n + 

5k,± = - ± -T7 

2 /i(cjk.- 



1)U - J k 



(50) 



(51) 



is the probability to create a particle (hole) excitation. 
Although the one-particle Green's function is a concept of 
importance in the context of quantum field theory, its use 
in the concrete experiments is limited by the impossibility 
to create the U(l) symmetry breaking interaction (|45[) . 
However, this function helps to interpret the nature of the 
excitation which is particle-like when one atom is added 
to the gas or hole-like when one atom is removed. 
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C. Sum rules 

Let us examine some sum rules satisfied by the suscep- 
tibility function. The compressibility sum rule is deduced 
from the Kramcrs-Kronig relation 



X "(k,a;) = x(k ) 0) 



(52) 



The /-sum rule generalizes the one obtained for a Bosc 
gas in continuum (37j 

h2 SZ 2^" (k '" ) = 
-^2 / d 3 k'cos(fc^)(a^ak')4Jsin 2 (fc Q /2) . (53) 



a=l 



In order to recover the continuum limit, we have to intro- 
duce the lattice constant a by means of the replacement 
k — > ak. In the limit of small a, we get [HI 



r k 2 



(54) 



where m = h 2 /(2Ja 2 ) corresponds to the effective mass. 
The application of Eq. (|52fl and Eq. (|5Uf in the Gutzwiller 
approximation leads to 



dw „, k^o 
— X (k,w) = -k 

1TUJ 



du 
2^ 



wx"(k,w) 



(0) 



Using the result (14TJ1) . we obtain 

Xk,A K 



E 



2 " 



(55) 
(56) 

(57) 
(58) 



A third sum rule concerns the static structure fac- 
tor S'(k). Using the fluctuations dissipation theorem, 
the dynamic structure factor is expressed as S(k, oS) = 
x"(k, ui)/it at zero temperature so that [l8j : 



/*°° A 

\ — X "(k,uj) = S(k) = (5p k 6p- k ) 

Jo 7T 



(59) 



where <5/3 k = X)i^i e-<k 'V id/2 - This third 

sum rule 

is not fulfilled in the Gutzwiller approximation because 
the correlation function is equal to the particle num- 
ber fluctuations (<5/3k<5/5-k) = (S 2 n) and thus has no 
k-dependence. However, as in the case of a Bose gas 
in continuum, this sum rule allows to deduce the static 
structure factor. We find indeed 



S(k) = ]>>,A k =°fc s |k 



(60) 



This last result shows an interesting feature of the sum 
rule approach. Starting from the lowest-order Gutzwiller 



approach that docs not contain any correlation, the 
two-point correlation function is determined as a next 
order contribution. Similarly, starting from the time- 
dependent DGPE, a similar procedure has been success- 
fully used to recover the static structure predicted from 
the Bogoliubov approach [TH, |60| • 



D. Spectra measurement 

The observation of the excitation branches in the SF 
phase can be made through the measurement of the total 
momentum P. After an adequate time of flight t, the 
momentum is given by [21 1 



P = J d 3 x^n(x) = hj d 3 kk|w(k)| 2 G(k) , (61) 

where w(p) is the Fourier transform of the Wannicr func- 
tion and 



G(k) = Y / e lW ' ) (dla l , 



(62) 



For small momentum, we can assume w(k) ~ w(0). Cal- 
culations up to the second order in the potential allow to 
deduce 



^ = -2kK0)| 2 |^| 2 Im X (k, W ) 



(63) 



= 27rk|w(0)-^| 2 ]T ±Xk,A<^ T Wk,A) ■ 



VI. CREATION OF SOUND WAVES 

Experimentally sound waves in a trapped Bose- 
Einstcin condensate were created turning on and off 
a perturbation potential in the center of the atomic 
cloud [6l|, [62| . The same can be done in optical lat- 
tices and numerical simulations of this kind of experi- 
ment were performed in Ref. [20[ deep in the SF phase 
making use of the DGPE and in Rcf. [14J for ID-systems 
using DMRG method. In this section, we do the same 
simulations but using the dynamical Gutzwiller ansatz. 
Our aim is to compare the results with that obtained by 
the other methods and to extract from the simulations 
the values of the sound velocity compatible with that 
calculated in Sec. IIVI 

We are interested in the solutions for d-dimcnsional 
lattices which have a position dependence only in one 
chosen spatial dimension. Then in Eq. © "0i±e„ = ipi±i, 
if a is the chosen dimension, otherwise ^i±e a = ipi ■ Here 
/ is the site index along the chosen dimension. Initially 
the atoms are prepared in the ground state of the external 
potential 



ei = e cxp [-(I - l ) 2 / 



1 I 21 

'w 



l = (£ + l)/2, (64) 



where £q and w are the strength and the spatial width 
of the potential, respectively The external potential 
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([64]) creates a density perturbation and we calculate 
the corresponding ground state numerically propagating 
Eq. ([6]) in imaginary time [(33| with the initial conditions 

Qn(0) = , where are the coefficients for the ground 
state of the homogeneous lattice with the local chemical 
potential /x; = /x — £;. At t = 0, the external potential 
([64]) is switched off (ei = 0) and the ground state starts 
to evolve. We calculate numerically the evolution of the 
ground state in real time. Numerical calculations pre- 
sented in this section are performed for d = 3 and we 
used L = 200, N = 10, which is enough to avoid para- 
sitic effects due to the reflection from the boundaries and 
due to the cut-off in the occupation numbers. 



A. Deep in the SF phase 

First we do simulations deep in the SF phase. In this 
case there is no difference between (hi) and \ipi\ 2 - Time 
evolution of (hi) for two negative values of £o is shown 
in Fig. [13] Initially the distribution of the atoms have 
a maximum at the position of the potential (bright per- 
turbation). After switching off the potential the density 
perturbation splits up in two wave packets propagating 
symmetrically outward. At longer times, the form of the 
propagating maxima becomes irregular which signals the 
creation of shock waves [l4|, H3] ■ The irregularities be- 
come more pronounced for larger values of |eo|- An addi- 
tional dip arises at the fronts of the wave packets which 
might stem from the discreteness of the lattice [l4[. For 
positive values of £o (gray perturbation), the dynamics 
is pretty much the same except that the maxima are re- 
placed by minima and vice versa and the distribution of 
the atoms is more regular (Fig. [14]). 

If we increase the width of the external potential w, 
the atomic distribution becomes more regular (compare 
Figs. nUT a) and [T5]) . In Fig. [IS] one can clearly see that 
the form of the wave packets become very asymmetric 
during their propagation. 

Figs . \W\ and [171 show the time dependence of the global 
maximum and minimum of the atomic distribution (hi) 
for negative and positive values of eq. When the external 
potential is switched off, the amplitude of the density 
perturbation goes down and after some finite time seen 
as the first plateau in Fig. [15] two separate wave packets 
are formed which propagate in opposite directions. Their 
amplitude decreases monotonically in time for negative £o 
(Fig.Qi)]) and for small enough positive £o (Fig. [T7J. For 
larger positive £o, the amplitude of the minima increases 
before it starts to decrease. Due to the discreteness of 
the lattice the position of the propagating cxtremuma is 
a step- function of time (Fig. IT5|) and the amplitude of the 
density perturbation shows up oscillations (Figs. [121 [T7j) 
which arc stronger for larger values of |eo|. 

Before the system enters the shock wave regime, there 
are always pronounced global maxima (minima) in the 
case of negative (positive) Eq and the propagation veloc- 
ity of the sound wave packets can be identified with the 
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FIG. 13. Time evolution of the mean occupation numbers 
(hi) deep in the SF phase after switching off the potential 
with w = 1, eo/U = —0.6 (a), —0.1 (b). The parameters are 
(h) — 1, 2dJ/U = 6. The curves show the spatial dependences 
at the dimensionless time r = tU/h = (0), 2 (1), 4 (2), 6 (3), 
8 (4), 10 (5). They are shifted by 0.4 (a) and 0.1 (b) in the 
vertical direction with respect to the previous one. 

velocity of the global extremum. Numerical values of the 
propagation velocity c s in our simulations are determined 
with the aid of a linear fit for the location of the global 
extremum as a function of time (straight lines in Fig. I18p . 
Its dependence on the amplitude of the external poten- 
tial is shown in Fig. 1191 Propagation velocity decreases 
monotonically with eq. This can be understood looking 
at the behavior of the function c®(n' = /j, — Sq) under vari- 
ation of £o at fixed 2dJ/U. For large values of 2dJ/U, 
it is also a decreasing function of £ (see Fig. [S] and the 
dashed line in Fig. [Hi]) , which is quite close to the data 
of our numerical simulations. 

In order to extract the value of the sound velocity from 
the numerical data, we have to extrapolate to £o = 0. 
This is done making a quadratic fit to the data points 
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FIG. 14. Time evolution of the mean occupation numbers 
(hi) deep in the SF phase after switching off the potential with 
w — 1, eo/U — 0.6. The parameters are (n) = 1, 2d J/U = 6. 
The curves show the spatial dependences at the dimensionless 
time r = tU/h = (0), 2 (1), 4 (2), 6 (3), 8 (4), 10 (5). They 
are shifted by 0.3 in the vertical direction with respect to the 
previous one. 
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FIG. 15. Time evolution of the mean occupation numbers 
{hi} deep in the SF phase after switching off the potential 
with w = 5, eo/U — —0.6. The parameters are (n) = 1, 
2dJ/U = 6. The curves show the spatial dependences at the 
dimensionless time r = tU/h = (0), 2 (1), 4 (2), 6 (3), 8 (4), 
10 (5). They are shifted by 0.3 in the vertical direction with 
respect to the previous one. 



which is justified by the fact that, near £o ~ 0, c° s (n' = 
(j, — eq) can be decomposed into a series in powers of e . 
In the example shown in Fig. [19l the extrapolated value 
of the propagation velocity is 1.433 which is a bit higher 
than Cg = 1.372 predicted by the linear response theory. 
We have done the same calculations for different values of 
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FIG. 16. Time evolution of the largest mean occupation 
number (hi)max after switching off the potential with w = 5, 
eo/U = -0.1 (1), -0.2 (2), -0.3 (3), -0.4 (4), -0.5 (5), 
—0.6 (6). The parameters are (n) = 1, 2dJ/U = 6. 
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FIG. 17. Time evolution of the smallest mean occupation 
number {hi) m in after switching off the potential with w = 5, 
eo/U = 0.1 (1), 0.2 (2), 0.3 (3), 0.4 (4), 0.5 (5), 0.6 (6). The 
parameters are (n) = 1, 2dJ/U — 6. 



w and found that the propagation velocity becomes more 
close to c° for larger values of w (see Fig. |20| . Therefore, 
the deviation from Eq. (|28| is due to the contribution of 
excitations with finite wavelengths. 



B. Near the boundary of the SF-MI transition 

Still in the SF phase but near the boundary of the 
SF-MI transition, \ipi\ 2 < (hi) (Fig [21]). Numerical sim- 
ulations in this regime show that, in order to excite only 



13 








2 4 6 8 10 

T 

FIG. 18. Location of the largest mean occupation number 
{ni)max after switching off the potential with w — 5, eo/U = 
-0.1 (1), -0.6 (6). The parameters are (h) = 1, 2dJ/U = 6. 
Here we use the same labels for the curves as in Fig. 1161 
Straight lines are linear fits to the numerical data. 
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FIG. 19. Dependence of the propagation velocity on the 
strength of the external potential. The parameters are (h) = 
1, 2dJ/U = 6, w = 5. The dots are the results of numerical 
calculations and the solid line is a fit by quadratic polynomial. 
The dashed line shows c°(// = fi — eo) as a function of Eo with 
u fixed by the values of (n) and 2dJ/U. 



the lowest mode, much less values of eo are required. This 
is consistent with the fact that the gap between the first 
and second excitation modes is very small near the MI 
lobe. 

Fig. [221 shows the dependence of the propagation ve- 
locity on Sq near the tip of the MI lobe with no = 1. It 
is quite different compared to the behavior deep in the 
SF regime where the propagation velocity monotonically 
decreases (Fig.QlJ]). Near the tip of the MI lobe, the prop- 
agation velocity has a maximum around eo ~ which is 
qualitatively similar to the behavior of c°(/i' = u — eo). 
However, the discrepancy between the propagation ve- 
locity and c°(/x' = /i — eo) is large even for small |eo| due 
to the significant contribution of nonlinear effects. 

We have determined the propagation velocity in the 
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FIG. 20. Dependence of the propagation velocity on the 
width of the external potential. The parameters are (h) — 1, 
2dJ/U = 6. 

limit eo — > making again a quadratic fit to the numer- 
ical data presented in Fig. [221 This procedure gives us 
the value 0.199, while from Eq. ([28]) we get c° = 0.201. 

C. MI phase 

If the parameters of the external potential e; are cho- 
sen such that the values of hi are always within the MI 
phase, the density is not perturbed and, therefore, there 
will be no time dynamics when the potential is switched 
off. In Fig. [231 we show an example of a broader po- 
tential, where local SF regions appear near the center of 
the potential. When the potential is switched off, the lo- 
cal inhomogeneities just spread without creation of any 
propagating wave packets. This is consistent with the 
fact that the sound waves do not exist in the MI phase. 

VII. CONCLUSION 

We have studied collective excitations of interact- 
ing bosons in a lattice at zero temperature within the 
framework of the gapless and conserving time-dependent 
Gutzwiller approximation. The excitation modes are 
calculated within the framework of the linear-response 
theory considering small perturbation of the many-body 
ground state. We demonstrated that the lowest-energy 
excitation of the SF has a phonon-like dispersion rela- 
tion and derived an analytical expression for the sound 
velocity in terms of compressibility and the condensate 
density which coincides with the hydrodynamic relation. 

We have studied the response of the lattice Bose gas 
in the Bragg scattering process which provides an ex- 
perimental tool to observe the excitations. It is demon- 
strated that the susceptibility function satisfies the f-sum 
rule in the whole parameter region. Calculations of the 
transition amplitudes show that within the Gutzwiller 
approximation the MI does not respond to the perturba- 
tion caused by the Bragg potential. In the SF phase, we 
show that only three lowest excitation branches have sig- 
nificant transition amplitudes, the others being too small 
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FIG. 21. Time evolution of the mean occupation numbers 
(hi) (a) and the mean number of condensed atoms \ipt\ 2 (b) 
near the tip of the Mi-lobe after switching off the potential 
with eo/U — —0.3, w — 5. The parameters are (n) = 1, 
2dJ/U — 0.172. The curves show the spatial dependences 
at the dimensionless time r = tU/h = (0), 24 (1), 48 (2), 
72 (3), 96 (4), 120 (5). They are shifted by 0.1 (a) and 0.2 (b) 
in the vertical direction with respect to the previous one. 
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FIG. 22. Dependence of the propagation velocity on the 
strength of the external potential. The parameters are (h) = 
1, 2dJ/U = 0.172, which is close to the tip of the MI lobe 
[2d{J/U)™ ax = 0.17157], w = 5. The dots are the results of 
numerical calculations and the solid line is a fit by quadratic 
polynomial. The dashed line shows c°(/i' = fi — eo) as a 
function of eo with /i fixed by the values of (n) and 2dJ/U. 



Appendix A: Derivation of Eq. (|28[) 

In order to work out the sound velocity c°, we con- 
sider the limit of small |k| and look for the lowest energy 
solution of Eq. f) 13[) as an expansion with respect to k: 



(Al) 



(7 k 


-,(0) 

= «k 






< 7 k 


-(0) 








- w k 




+ 4 2) 



The zeroth-order solution satisfies the equation 

,(°) 



a b W4°) 



(A2) 



, u k / v " " > \ L 'k / 
and is non-trivial only in the SF phase with the form 
dTiLOi 



,(°) - 



dfx 



J ° \ Jo) 



v (0) _ (0) (0) Q 



(0) 
kn ' 



(A3) 

The quantities of the first order are governed by the equa- 
tion 



to be observed in a real experiment. The absence of re- 
sponse in the MI phase is a limitation of the Gutzwiller 
approximation and can be corrected in the next leading 
order in the inverse of the coordination number which 
would take into account the possibility of a particle-hole 
pair creation (2f| H3] • 

Finally, we have performed simulations of the sound- 
wave propagation solving numerically the Gutzwiller 
equations. The calculations show that sound waves can 
be created only in the SF phase and the corresponding 
velocity is in a good agreement with the results of the 
linear-response theory. 



Taking into account the identity 

E(^'+*r')^ = ( n+ ^U 0) ' (as) 

n' r / 

the first-order solution can be written as 

«S=«£!=^ ) %-. (A6) 
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Multiplying Eq. (|A7[) by the vector (u^ , — v^^) from 
the left side and taking into account that, 



(o) (o) 

— V 



kn ' dfx dn 



8(h) 

— K > 



(A8) 

(A9) 
for 



where k is the compressibility, we arrive to Eq. 
the sound velocity. 



Appendix B: Bogoliubov's dispersion relation 

We look for the solution of Eq. ([TU| for arbitrary k in 
the form 



AO) 



(Bl) 



where u^ 1 ^ and v^ 1 ' are given by Eqs. (|A3|) . (|A6|) with 



(0,1 



X 1 ) 



being replaced by un- 
plugging (|B1[) into Eq. ([l"3]) and multiplying the re- 



sulting equations by vectors u k , u^ 1 ' 
homogeneous equations for a k and fe k : 



we obtain linear 



E 
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FIG. 23. Time evolution of the mean occupation num- 
bers (hi) (a) and the mean number of condensed atoms 
\ipi\ 2 (b) in the MI phase after switching off the potential 
with Eo/U — 1, w — 5. The parameters are iijU = 1.2 and 
2dJ/U = 0.07, which is a bit smaller than the critical value 
2d(J/U) c — 0.073. The curves show the spatial dependences 
at the dimensionless time r = tU/h = (0), 20 (1), 40 (2), 
60 (3), 80 (4), 100 (5). They are shifted by 0.5 in the vertical 
direction with respect to the previous one. 



We substitute all these results in the equation for the 
quantities of the second order in k 




(A7) 



Ht "k huJkK/2 



2 + \~djT 



<~k 



E4 2(^r'-sr')4°^ = ^ 



(0)' 



lead us to the result 



\ 




(B2) 



£«£(41 n ' (B3) 



(B4) 



(B5) 



In the limit of small U /J, (h) = ip^ and we recover the 
well-known Bogoliubov's dispersion relation. For arbi- 
trary values oil// J, Eq. (|B5|) better describes the lowest 
excitation branch of the SF than the standard Bogoli- 
ubov's dispersion relation but gives higher values of en- 
ergies compared to the exact numerical data. 
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Appendix C: Energy gaps 

We first rewrite Eq. in the form 

where the operators 

^4k = —Ja^°\a + a)) + -^n(n — 1) — fin — huo 
B k = - J k a|sW) (s<°> |a + at| S (°))( S (°) |a+ , 

with the ground state \s^) defined by Eq. ([3|), act on 
the kets 

\uk) = ^u kn \n) , \vu) = y^fknN . (C2) 

n n 

The zeroth order exact solution in the small parameter 



U I J is given by 

i«S> - 6 ^ (0) ) i A ) . 

|4a) = , (C3) 

where -D(a) = exp(ad^ — a* a) is the displacement op- 
erator. This can be shown using the property of the 
displacement operator 

D(a)ab{a) = a + a . (C4) 

Finally, first order perturbation theory in small U/J and 
the relation 

<«Sl3k|«i3> = ( C5 ) 
valid for A > 2 allow to establish that 

^£1 = (ui°l\^\u£l) . (C6) 

Explicit calculation of the matrix element leads to the 
result |27|) for the gaps in the excitation spectrum. 
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